function [Yp,fai] = ccm(L,Y)
n = size(L,1);

cls = unique(Y);
c = size(cls,1);

disp('Start optimization...');
tic
options = optimoptions('quadprog','Algorithm','trust-region-reflective');

fai = zeros(n,c);

for i = 1:c
    disp(['Training class: ',int2str(i)]);
    ipos = find(Y==cls(i));
    lpos = length(ipos);
    ineg = find(Y~=cls(i));
    lneg = length(ineg);
    Yl = -inf*ones(n,1);
    Yu = inf*ones(n,1);
    Yl(ipos) = ones(lpos,1);
    Yu(ineg) = -ones(lneg,1);
    x0 = zeros(n,1);
    x0(ipos) = ones(lpos,1);
    x0(ineg) = -ones(lneg,1);
    fai(:,i) = quadprog(L,[],[],[],[],[],Yl,Yu,x0,options);
end

[~,Yp] = max(fai,[],2);
Yp = cls(Yp);
end